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Abstract: We analyse the associated production of Higgs and Z boson via heavy-quark 
loops at the LHC in the Standard Model and beyond. We first review the main features of 
the Born 2 —>■ 2 production, and in particular discuss the high-energy behaviour, angular 
distributions and Z boson polarisation. We then consider the effects of extra QCD radia¬ 
tion as described by the 2 —3 loop matrix elements, and find that they dominate at high 
Higgs transverse momentum. We show how merged samples of 0- and 1-jet multiplicities, 
matched to a parton shower can provide a reliable description of differential distributions in 
ZH production. In addition to the Standard Model study, results in a generic two-Higgs- 
doublet-model are obtained and presented for a set of representative and experimentally 
viable benchmarks for Zh °, ZH° and ZA° production. We observe that various inter¬ 
esting features appear either due to the resonant enhancement of the cross-section or to 
interference patterns between resonant and non-resonant contributions. 
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1. Introduction 

Run I of the LHC has brought the discovery [1, 2] of a scalar particle, so far consistent 
with the Higgs boson predicted by the Standard Model (SM) [3], and it has given the first 
evidence of the Brout-Englert-Higgs mechanism [4, 5] in particle physics. The increased 
energy and luminosity that will be achieved in Run II at the LHC will allow us to pin down 
the properties and in particular the strength and form of the interactions of such a boson 
with all other SM particles. To this aim a vast campaign of measurements of rates and 
distributions in various production and decay channels is being planned. 

Among such processes is the associated production of a Higgs boson together with a 
vector boson V, either a or a Z, also known as Higgs-strahlung, i.e., at the leading order 
in QCD, the Drell-Yan production of an off-shell vector boson qq —>• V* with its subsequent 
decay V* —>• VH. While suppressed in the SM with respect to the leading gluon-gluon 
and vector boson fusion channels, VH production is of phenomenological interest mostly 
because the presence of the vector boson (and possibly of leptons coming from its decay) 
in the final state can help to access the large yet challenging H —» bb decay mode. For 
instance, Higgs-strahlung has been the dominant Higgs search mode at the Tevatron [6]. At 
the LHC, the ATLAS [7] and the CMS [8] collaborations have investigated VH production, 
with the Higgs boson decaying to a b —quark pair, both reporting small excesses above the 
background only hypothesis. Searches for Higgs decaying to W + W~ [9,10] and to invisible 
states [11,12] have also been performed by both ATLAS and CMS. 

On the theory side, predictions for ZH production are known at NNLO in QCD and 
at NLO electroweak in EW theory. The NNLO QCD cross section includes the Drell- 
Yan type terms of (D(g 4 a 1 2 s ) first computed in [13,14]. In addition to Higgs-strahlung, it 
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has been noted that contributions from quark-anti-quark initiated diagrams where the 
Higgs is emitted from a top quark loop arise at the same order. These diagrams inter¬ 
fere with the LO and NLO Drell-Yan amplitudes and have been computed in [15], where 
they were found to contribute to the inclusive NNLO cross section at the percent level. 
Implementations of the NNLO QCD calculations are publicly available in vh@nnlo [16] 
and HVNNLO [17]. Fully differential NLO QCD and EW results can be obtained with 
the program HAWK [18,19], while event generation accurate at NLO in QCD (inclusively 
and for higher jet-multiplicitites), can be nowadays obtained (automatically or semiauto- 
matically) in several frameworks, i.e., MadGraph_aMC@NLO [20] / POWHEG [21] + 
Pythia 8 [22] /HERWIG++ [23] and Sherpa [24], 

At NNLO, a purely virtual gluon fusion contribution emerges, through the gg —>• ZH 
amplitude squared, which at the LHC can be enhanced by the large gluon-gluon luminosity 
at small Bjorken x. Its contribution to the total cross section has been known for a long 
time [25,26] and it has been included in the implementations of the NNLO calculations [16, 
17]. The gluon fusion component is separately gauge invariant, IR and UV finite and 
accounts for about 10% of the total NNLO cross section at 14 TeV. Being essentially a 
leading-order contribution, gg —> ZH introduces a rather strong scale dependence to the 
NNLO result, which in turn is known quite precisely. In order to reduce the associated 
theoretical uncertainty, recently, NLO corrections for the gluon fusion contribution have 
been estimated by computing them in the infinite top-quark mass limit [27]. The NLO 
corrections to this process, O(af), while formally part of the N 3 LO ZH cross section, 
are expected to be large, similarly to other gluon fusion processes such as Higgs single 
or pair production. The computation of the approximate NLO corrections in the infinite 
top mass limit has confirmed this expectation. The NLO computation in the infinite 
top mass limit reduces the scale uncertainty by a factor of two, yet the size of the finite 
top-quark mass effects remains unknown: the exact NLO result requires two-loop multi¬ 
scale amplitudes whose analytic form is beyond the current advances in the multi-loop 
technology. In an effort to further reduce the theoretical uncertainties in this process, a soft 
gluon resummation for the gluon-gluon contribution has been performed in [28] promoting 
the previous results to NLO+NLL accuracy. We should note here that in contrast with 
single Higgs production, where the infinite top-quark mass limit provides a good description 
of the process, and allows the computation of higher order corrections, here, similarly to 
(yet with even less control than) Higgs pair production, the much higher scales involved 
make the effective field-theory approach unreliable, especially so at the differential level. 

In addition to the SM production mechanism and characteristics, interesting features 
can be expected from Higgs production in association with a Z boson in beyond the SM 
scenarios. The Two-Higgs-Doublet-Model (2HDM) is an attractive framework in which 
Higgs-strahlung can lead to interesting features. First the range of channels is richer: in 
addition to the production of the light (125 GeV) Higgs boson in association with a Z 
boson ( Zh °), Z associated production of the heavy scalar ( ZH °) and pseudoscalar boson 
(ZA°) are also possible [29]. Experimental searches are already underway to look for 
signals of these processes, especially in the case where the cross-sections can be enhanced 
by the resonant production of an intermediate scalar (H° or A 0 ) with subsequent decays 
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into Z and a lighter scalar. In particular, CMS has searched for signals of the decay of the 
pseudoscalar A 0 into a Zh° pair [30] and that of the heavy scalar H° into a ZA° pair [31] 
and the results have been used to set constraints on the 2HDM parameter space. 

So far considerable effort has been devoted to provide accurate total rates in this 
channel for both the SM and the 2HDM, but accuracy and precision in the differential dis¬ 
tributions is also of vital importance. This need becomes more important for experimental 
analyses which make use of exclusive observables, in order to tame the typically very large 
QCD backgrounds. Moving in that direction, it has been noted in the literature [29, 32] 
that the gluon induced component can play an important role. The gluon fusion Higgs 
Pt distribution peaks at higher values than the corresponding Drell-Yan one, and there¬ 
fore its relative contribution becomes more important in boosted Higgs searches, which 
are preferred experimentally to reduce the backgrounds. The prospects of such searches 
have improved recently due to progress in jet substructure techniques, after the seminal 
suggestion in [33]. 

The aim of this work is to contribute to the understanding of gluon induced ZH 
production and to improve the predictions for the differential distributions. We consider 
the 2 —s- 2 and 2 —> 3 matrix elements entering the gluon fusion contribution to Z Higgs- 
strahlung. We first review the main features of the 2 —>■ 2 ones and then examine the 
importance of the 2 —> 3 contributions. Given the lack of an exact and fully differential 
NLO computation for this process, we provide a better description of the kinematics for 
this component by combining the 2 —>• 2 and 2 —*• 3 matrix elements in a merged sample, 
matched to a Parton Shower (PS). This provides a fully exclusive control at the hadron 
level. A similar approach has been followed for other loop induced processes, such as single 
Higgs production [34] and Higgs pair production [35,36]. In general, this method provides 
a better description of the kinematics, yet as the formal accuracy for total rates remains at 
LO, it is often combined with a normalisation obtained from higher-order computations, 
when available. 

This merging-matching approach makes use of the fact that while tree level fixed- 
order amplitudes describe reliably the region of hard and well separated jets, the parton 
shower provides a better description of the soft and collinear regions. Combining the 
two requires of course a consistent treatment to avoid double-counting, which is achieved 
by various merging algorithms. Methods that are widely used for tree level merging are 
CKKW [37,38], CKKW-L [39,40] (and their later improvements [41,42]), and MLM [43]. 
More recently new methods have been developed to perform the merging at NLO, see for 
example FxFx [44] and UNLOPS [45], yet not directly applicable to 2 —>• 2 loop-induced 
processes at the Born level yet, mainly due the absence of analytic results for the two-loop 
2 —>■ 2 matrix elements. 

In this work we study gluon induced Higgs-strahlung at the LHC, presenting the first 
merged-matched results for gg —>• Zq i>, with (j) being a generic scalar, by employing the 0 
and 1-jet matrix elements for the SM and the 2HDM. This paper is organised as follows. 
In Section 2, we discuss the process within the SM, first by reviewing the important fea¬ 
tures coming from the gg —>• ZH matrix elements. We also consider the behaviour of the 
2 —>■ 3 matrix elements, which we then combine with the 2 —>• 2 ones. We describe our 
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methodology, and present results both at the parton level and after merging and matching 
to a parton shower. In Section 3, we explore the results of various 2HDM scenarios using 
the same calculation setup. We draw our conclusions in the final section. 

2. Gluon induced ZH production in the SM 

Representative Feynman diagrams contributing to the gg —» ZH process in the SM are 
shown in Fig. |l]. Massive fermions, t and b— quarks, run in the box, while all flavours run 
in the triangle. The contribution of the two light generations to the triangle vanishes as 
required by the anomaly cancellation. In practice, it is only the axial vector part of the 
heavy-quark-Z coupling that contributes to the amplitude. The amplitude for this process 
was first computed in [25,26]. 

g 

mmm? 

g 

Figure 1: Representative Feynman diagrams for ZH production in gluon fusion in the SM. 

In what follows, we will first review the main features of the 2 —>■ 2 process for gluon 
induced ZH production before discussing the implications of the 2 —>■ 3 one. A sample 
of the relevant diagrams contributing to ZHj is shown in Fig. ||. In addition to the gg 
initial state amplitudes, the qg and qq channels also open up, when an additional jet is 
allowed. The gg —>• ZHg amplitudes were used in [46] to calculate the gg part of the ZHj 
cross-section at the LHC for various jet transverse momentum cuts. In what follows, we 
will consider these along with the qg and qq diagrams to discuss the behaviour of the 2—^3 
amplitudes and subsequently to obtain a merged sample of 0 and 1-jet multiplicitities. 

2.1 Calculation setup 

In this work, we employ the MadGraph5_aMC@NLO framework [20]. The one-loop 
amplitudes squared for ZH and ZHj can be obtained with the help of MadLoop [47], 
which computes one-loop matrix elements using the OPP integrand-reduction method [48] 
(as implemented in CutTools [49]). A reweighting procedure is then employed to over¬ 
come the present limitations concerning event generation for loop-induced processes 1 . A 
reweighting method has been employed already for a series of processes within the Mad- 
Graph5_aMC@NLO framework [34,51,52] both at LO and NLO accuracy. This procedure 
involves generating events through the implementation of a tree-level effective field theory 
(EFT), in this case obtained by taking the limit of infinite top-quark mass with all other 
quarks being massless. In practice, a UFO model [53,54] including the effective theory in¬ 
teractions is imported in the simulation framework. After event generation, event weights 

1 Automated event generation for loop-induced processes is currently being finalised [50]. 
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Figure 2: Representative Feynman diagrams for gluon induced ZHj production in the SM. 

obtained from the tree-level EFT amplitudes are modified by the ratio of the full one-loop 
amplitude over the EFT ones, i.e., r = \M 2 Loop \/\M 2 EFT \, where \M^\ oop \ represents the 
numerical amplitude as obtained from MadLoop. In our case, reweighting proves to be 
efficient in terms of the computational speed, as the loop amplitudes have to be calculated 
for significantly fewer phase-space points than what would be needed to integrate them 
directly. Moreover the EFT leads to distributions that are in general harder in the tails, 
and therefore the EFT events populate regions that are later suppressed by the exact loop 
matrix elements, resulting to no significant degradation of the statistical uncertainty. 

2.2 Parton level results 

Before proceeding to the technical setup and presenting results of the merging-matching, 
we consider the salient aspects as observed at the parton level. The findings of this study 
will reveal some previously unnoticed features of gg —>• ZH and will act as a motivation to 
employ a merging-matching procedure in the following section. 

In our computation the heavy quark masses are set to: mt =173 GeV and m& =4.75 
GeV, while the Higgs mass to mn =125 GeV and the heavy quark Yukawas are given by 
y q /V2 = rriq/v. We note here that finite width effects in the propagators of the loops can be 
taken consistently into account within MadGraph5_aMC@NLO via the implementation 
of the complex mass scheme [55,56]. The effect of a non-zero top width is shown in Fig. |3], 
where the matrix element squared for gg —» ZH , for 90° scattering, is shown as a function 
of the invariant mass of the ZH system. The correction is more important at the tt 
threshold, where it reaches 20%. Finally, when integrated over all centre-of-mass energies 
and scattering angles, we find the top-quark width to modify the gg —» ZH cross-section 
by rs_/ 2% at 14TeV, an effect similar to that observed for single and double Higgs production 
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in [57] and [52], respectively. For the rest of the results presented in this work the width 
of the top quark is set to zero. 



300 400 500 600 700 800 900 1000 

s/S [GeV] 


Figure 3: Top width effect on the matrix element squared for gg —> ZH. Results for 1^=0 and 
1.5 GeV are shown along with their ratio. 

An interesting aspect of the gg —)• ZH matrix element is its angular dependence. While 
in Fig. H] we have fixed the scattering angle to 90°, in Fig. |||, we show the dependence of 
the amplitude squared on the centre-of-mass scattering angle, for various values of s. The 
matrix element starts with no angular dependence at low energies, but varies significantly 
with the angle at high energies. This angular dependence of the matrix element implies 
that at high energies, very forward or backward scattering is favoured over 90° scattering. 
This behaviour originates from the interplay between the triangle and the box diagrams, 
and their respective angular behaviour. As we will also discuss later, box and triangle 
interfere destructively, with the triangle contribution dominating at low energies. The 
cancellation becomes nearly exact at high energies, mostly leaving a remainder from the 
box contribution that is strongly dependent on the scattering angle. 

We now proceed to discuss results for the LHC. For these, parton distribution functions 
(PDFs) are evaluated using the MSTW2008LO set [58] and the central renormalisation 
and factorisation scales are set to the invariant mass of the ZH system: pP = pP R = pP F = 
mzH■ MadGraph5_aMC@NLO allows the automatic computation of the scale and PDF 
uncertainties by a reweighting procedure [59]. In our results, scale variations are obtained 
by varying the scales in the range of p° / 2 < pr : f < 2/x°. 

Table |] summarises results for the gg —» ZH cross section and the “loop-induced” ZHj 
contribution, originating from the square of the amplitudes shown in Fig. |j. For reference 
we also include the total pp —>• ZH cross section at NNLO obtained with VH@NNLO [16]. 
The NNLO cross section includes the gg —>• ZH result of the first row, for which excellent 
agreement has been found between our computation and the result of VH@NNLO. We note 
that the results shown in the second row of Table [fj for ZHj are obtained using the loop 
amplitudes shown in Fig. |j. These qg and qq amplitudes can interfere with the Drell-Yan 
type real emission amplitudes. This interference contribution to the cross section has been 
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Figure 4: Matrix element squared for gg —> ZH as a function of the centre-of-mass energy for 
different values of scattering angles (left) and as a function of the centre-of-mass scattering angle 
for various values of y/§ (right). 


computed in [15] and found to be at the per-mille level. In our computation we use these 
amplitudes squared, i.e., at O(a^). It is clear that at this order, other qg loop-induced 
contributions can enter squared, for example, the set of diagrams where the Z couples 
to a light quark and the H to a top-quark loop. We have not included these diagrams 
here, as we consider them of a different origin, but we have checked that their amplitude 
squared contribution to the cross section is small, below the femtobarn level, and therefore 
at least one order of magnitude smaller than those in Fig. ||. The interference of this type 
of diagrams with the Drell-Yan amplitude was computed in [15] and also found to be small. 

Given that the ZHj amplitude is divergent in the limit of a collinear or soft jet, we 
apply a cut on the px of the jet to obtain finite results. We have set this cut to 30 GeV 
in Table [I]. The 2—^3 contribution comes mainly from the gg initiated diagrams, with qg 
giving about 20% of the ZHj cross section. The ZHj contribution is not as suppressed as 
expected from the extra power of a s , leading to results comparable in size to the gg —>• ZH 
cross section. Of course these results are extremely sensitive to the chosen cut for the 
transverse momentum of the additional parton, as the cross section diverges in the IR 
limit. Such a problem would not arise in the case of a NLO computation matched to a 
parton shower, for example with the MC@NLO method [60], in which all divergences are 
regularised and cancelled for inclusive observables. 

The results in Table |] also demonstrate the problem of large scale uncertainties for the 
LO gg cross section, that contribute significantly to the total NNLO scale uncertainty. The 
problem persists also for the loop-induced ZHj contribution. A significant reduction of 
these intrinsic QCD uncertainties can only be achieved by a complete NLO computation, 
as discussed in [27]. 

In addition to the total cross-section results presented in Table [|, interesting observa¬ 
tions can be made by studying the differential distributions for the gluon fusion process. 
We start by presenting distributions for the invariant mass of the ZH system and the 
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Contribution [fb] 

Vi = 8 TeV 

yfs = 13 TeV 

Vi = 14 TeV 

gg -+ ZH 

pp —> ZHj (pip > 30 GeV) 
pp -+ ZH (NNLO) 

17 4 + 34 % 

1 ' -24% 

i 9 A +52% 
iz.u _ 32% 

007 + 2 . 2 % 

00 ' - 1 . 6 % 

co c +30% 
- 21 % 

49 0 +44% 

UJ.U _ 32% 

7Q7 +3.2% 

' y0 - 2 . 0 % 

70 n +29% 
- 21 % 

CO A +47% 
58.4 _ 31% 

QQC +3.2% 

88b _ 2 _ 3% 


Table 1: Cross sections (in fb) for ZH associated production at the LHC at •+ = 8,13 and 
14 TeV. The uncertainties (in percent) refer to scale variations. No cuts are applied to final state 
particles apart from the jet pr cut in the second row (jP T > 30 GeV) and no Higgs or Z branching 
ratios are included. The ZHj contribution shown here is obtained from the loop diagrams shown 
in Fig. 2, while the NNLO results are obtained with vh@NNLO [16]. 
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Figure 5: Invariant ZH mass and + distributions for ZH production at yfs = 14 TeV. The gluon 
fusion contribution is decomposed into the triangle and box contribution. For completeness we also 
plot the Drell-Yan type contributions at NLO obtained with MadGraph5_aMC@NLO. 


transverse momentum of the Higgs in Fig. ||| for the LHC at 14 TeV. These distributions 
have been shown elsewhere in the literature, for example in [32], yet we consider them 











































here again for completeness. In addition to the gluon fusion results, we also include the 
NLO Drell-Yan like distributions obtained automatically with MadGraph5_aMC@NLO. 
We note here that differential NNLO results for the Drell-Yan like contribution can be 
provided by the code HVNNLO. As is evident from [17], the NNLO computation leads 
to an overall 20% decrease of the Drell-Yan component but not to any significant shape 
difference compared to the corresponding NLO one. 

The first observation regards the clear presence of the 2 mt threshold, at which the 
gluon fusion amplitude acquires an absorptive part, related to the on-shell gg —>• tt , 
tt —>• ZH scattering, leading to a characteristic rise in the invariant mass distribution. It is 
evident from Fig. || that the gluon fusion component leads to distributions of fundamentally 
different shape from the Drell-Yan ones and therefore it should be considered in all relevant 
studies, in particular in the boosted region of p^ > 100 GeV, where its relative importance 
increases. 

In the plots, we decompose the gluon fusion result into the triangle and box compo¬ 
nents. The two interfere destructively over the whole range of centre-of-mass energies, with 
the cancellation between the two being nearly exact at high energies. Such a cancellation 
is due to unitarity: while each of the two diagrams grows with energy the cancellation 
leads to a well-behaved amplitude at high energies. We stress here that this behaviour of 
the amplitude is not present in the infinite top mass limit. In this limit, the amplitude 
for the box diagram vanishes and therefore only the triangle contributes to the amplitude, 
giving a rather bad approximation of the one-loop amplitude at high energies. In addition, 
we note that this is a process highly sensitive to the relative phase between the HZZ and 
Htt couplings. To demonstrate this, in Fig. |] we show the result obtained by changing the 
relative sign between the top Yukawa and the HZZ coupling. In pair with other processes 
where such unitarity cancellations take place, such as H —>• 77 or pp —>• tHj [61-63], flip¬ 
ping the sign results in an increase in the gluon fusion induced contribution by a factor of 
five, and much harder distributions as the interference between triangle and box becomes 
constructive, see Fig. [5]. We conclude that, given the difference in the shape as well as the 
size of the cross section above 2 mt, the ZH invariant mass or transverse momentum of the 
Higgs or the Z could also be used to bound the relative phase between the Higgs couplings 
to fermions and to vector bosons. 

The difference in the px shape between the Drell-Yan and gluon fusion production 
persists also in the distribution of the lepton px coming from the Z decay. Besides, an¬ 
other interesting aspect of the gluon fusion process is that it leads to different angular 
distributions for the resulting leptons compared to the Drell-Yan component [64], This is 
evident from studying the normalised distributions of the angle 9*, defined as the angle of 
the lepton between the lepton and Z direction in the Z rest frame is shown in Fig. ||. In the 
plot, we use the NLO Drell-Yan result, and plot the distributions with and without a cut of 
100 GeV on the px of the Z. The shape of the distribution without any cut, is significantly 
different, with the tree-level ZH giving a flat distribution while the gluon-fusion one peaks 
at 90°. The shape becomes similar for pj, > 100 GeV, while a 200 GeV cut (not shown 
here) completely eliminates the difference. This behaviour is related to the polarisation of 
the Z which differs between the two production modes. This can be quantified by examin- 
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Figure 6: Normalised distributions of the lepton angle 0* for gluon fusion and Drell-Yan like ZH 
production at yfs = 14 TeV. Results for an imposed cut of 100 GeV on the pp of the Z are shown 
in dashed lines. 


Process 

/o (%) 

h (%) 

}r (%) 

gg ->• ZH 

82.2 

8.9 

8.9 

gg —> ZH , pj, > 100 GeV 

86.3 

6.9 

6.8 

qq —>• ZH 

35.6 

32.4 

32.0 

qq —> ZH , p^ > 100 GeV 

62.6 

18.8 

18.6 


Table 2: Polarisation fractions for the gluon fusion and NLO Drell-Yan production mode of ZH 
at 14 TeV with and without a cut on the Z pp- 


ing the relevant polarisation fractions in Table |2], as these are defined in [65]. The fact that 
the Z in gg —>• ZH is predominantly longitudinal leads to the central peak, while the small 
difference between fp and fp leads to a very mild asymmetry for qq. Setting a 100 GeV 
cut on the Z pp changes these values in agreement with the equivalence theorem, he., by 
increasing the longitudinal polarisation fraction. For completeness, we also mention here 
that the main background for this process, Z + 6-jets leads to predominantly left-handed 
Z bosons [65,66] and therefore to different angular distributions, that could be potentially 
used as an additional discriminating handle to distinguish signal and background. 

Further to the gg —> ZH results that we have discussed above, interesting conclusions 
can be drawn by studying the loop-induced ZHj distributions. We have seen in Table 
[l] that these contributions are not negligible and their relative importance increases with 
the centre of mass energy. A complete NLO computation for gg —» ZH would be fully 
inclusive in these contributions but as such a computation is not available yet, we aim to 
draw some conclusions by studying them independently. For such a study a minimum cut 
has to be set on the transverse momentum of the additional jet to avoid the divergent soft 
and collinear limit. We compare the distributions of the invariant mass of ZH and the p^ 
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to those from 2—^2 amplitudes by varying the px cut set on the additional jet in Fig. |7[ 




0 200 400 600 800 1000 

Pfi [GeV] 


Figure 7: Invariant ZH mass and ptf distributions for loop-induced ZH and ZHj production at 
yfs = 14 TeV. The results for the ZHj distributions are shown for various jet px cuts: 30, 50 and 
100 GeV and again concern only the loop diagrams of Fig. 2. 


The ZH invariant mass distribution shows that for all values of the jet px cut the 
bulk of the cross-section remains close to the 2m* threshold. The characteristic threshold 
behaviour due to the absorptive part of the amplitude remains visible for all cuts. At high 
invariant masses, we find that the amplitudes with an extra parton fall more slowly and 
overtake the 2 —2 process. 

In contrast to the invariant mass distribution where no extreme modification of the 
shape takes place, the p^ distribution is very much affected. First, we note that the 
threshold corresponding to mzH ~ 2 nit is now not visible in the px distributions for ZHj. 
The second and more striking observation is that above 300 GeV the 2 —>■ 3 process leads 
to a much harder px spectrum compared to the 2—^2 one. Moreover, for pif values above 
400 GeV all three distributions for ZHj coincide. The fact that in this region the result 
is insensitive to the jet px cut implies that hard jet emissions are dominating. This occurs 
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because by allowing the emission of a parton new kinematic configurations open up. In 
this high pt region, the kinematic configuration in which a soft jet is emitted and the Z 
and H basically recoil against each other is not the most favourable one. Instead, the 
configuration in which a hard jet recoils against the H, with the Z remaining rather soft 
becomes the preferred one. We have explicitly confirmed this behaviour by setting a high 
cut on the px of the Higgs, and studying the corresponding jet and Z transverse momentum 
distributions. A clear preference for the configurations where the jet is hard and the Z 
is rather soft is found when sufficiently far from the IR divergence. The behaviour of the 
2 —>■ 3 amplitudes at high px can be traced back to the presence of t —channel gluon diagram 
such as the gg —>• ZHg one shown in the top right of Fig. which becomes dominant in 
this region. The same behaviour is displayed by the qg —>• ZHq contributions, when these 
are considered separately, as they include diagrams of the same type as shown in the second 
row of Fig. ||. 

In conclusion, we have found that, especially for the transverse momentum distribution 
of the Higgs, the emission of an additional jet can dramatically modify the shape, due 
to new allowed kinematic configurations. This effect might prove important in studies 
involving highly boosted Higgs as discussed for example in [32]. The 2 —>■ 3 matrix elements 
are important and therefore need to be taken into account for accurate simulations. To 
combine the two in a consistent way and therefore provide a more realistic picture of the 
differential distributions, we will resort to merging and matching to a parton shower. In 
the following section we will discuss how this method allows us to provide more accurate 
predictions for the distributions. 

2.3 Merging different jet multiplicities: setup 

Given the lack of a complete NLO computation and the relevance of the 2 —>• 3 matrix 
elements, the best available procedure to accurately predict the distribution shapes is to 
employ the Matrix-Element-Parton Shower (ME+PS) procedure. ME+PS schemes allow 
the consistent combination of matrix elements with different jet multiplicities via their 
matching to a parton shower. In our study we employ the fcr-MLM scheme as implemented 
in MadGraph5_aMC@NLO. Merged samples are then passed to Pythia 8 [22,67] for 
PS. 

The implementation of MLM merging in MadGraph5_aMC@NLO/Pythia 8 comes 
in two variants: the traditional fcr-MLM and the shower-fc^ schemes. The two give com¬ 
parable results as discussed in [68]. In this study we will employ the shower-fc^ scheme. 
While this scheme has been used for phenomenological studies with Pythia 6 in the 
past, see for example [69], we employ the most recent implementation of the scheme in 
MadGraph5_aMC@NLO combined with Pythia 8. For a detailed description of this 
approach one can refer to [68], while here we only mention its main points. 

In the shower -kx scheme, matrix element events are generated with a minimum separa¬ 
tion pr min) between parton and the initial state ( iB ), and Q cn t between final-state partons 
(ij), defined by the measure: 

dw = VXi > Pr min , 4 = m HPTi’PT^ AR h > ^ut> (2- 1 ) 
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where A Rfj = 2[cosh(?yj— r]j)~cos(fa—4>j)\ and pr { , Vi and fa are the transverse momentum, 
pseudorapidity and azimuthal angle of particle i. Short distance (parton level) events 
are then passed to Pythia 8 which evolves them down using the p^-ordered shower. In 
practice, for each event Pythia 8 records the scale of the hardest shower emission: Q^ardest- 
This scale is then used to accept or reject the event as follows: for the low multiplicity 
events, the event is rejected if ^hardest > <3cut> while for the highest multiplicity the event 
is rejected if Qhardest > Q softest > whh Q^f est being the scale of the softest matrix-element 
parton in the event. The value of Q cnt is selected on a process-by-process basis to ensure 
that there is a smooth transition between the ME and PS regimes. In practice, this is 
assessed by examining the differential jet rate distributions which show if the transition is 
indeed smooth. 

2.4 Merged sample results for ZH in gluon fusion 

Using the setup described in the previous subsection for the merging and matching, we 
present in this section our merged results for various observables. In our simulations we 
keep the H and Z stable. For the merging performed here, the shower-A:^ scheme is used 
with Qcut = PT min = 30 GeV. We have checked that this choice leads to smooth differential 
jet rate distributions, and therefore a smooth transition between the ME and PS regimes. 

We start by presenting the results for the invariant mass of the ZH system and the px 
of the Higgs in Fig. ||[ while p^ H and distributions are shown in Fig. |9|. A comparison 
is made between the gg — » ZH sample showered with Pythia 8 and the merged 0 and 1- 
jet matched sample, presented in combination with the uncertainties associated with scale 
choices for both the factorisation/renormalisation scale of the hard process and the shower 
starting scale. We set the central value for the renormalisation and factorisation scales to 
mzH , as for the parton-level results. The shower starting scale in Pythia 8 can be set 
to either the kinematical limit (px = ^r), corresponding to what we refer to as “power” - 
shower or the factorisation scale of each event (mzH in our case), i.e., “wimpy”-shower. 
Pythia 8 allows us, for the “wimpy”-shower case, to modify the shower starting scale in 
the range of 0.5 pp < Qps < 2 pp . This gives us the possibility to systematically study 
the dependence of the results on the choice of the shower scale for both the merged and 
gg —> ZH- only samples, as shown by the blue bands in the plots. To study the systematic 
uncertainties due to renormalisation and factorisation scale variations, we vary the scales 
between 0.5/i° < pr.f < 2/y°, with /j° = mzH- This variation is shown by the yellow bands 
in the plots (with the central prediction being the “power”-shower result). In the results 
shown for mzH and Px hr Fig. we also include the parton-level results for comparison 
purposes. 

First, we notice that not all distributions are sensitive to the procedure of merging¬ 
matching. In particular, the invariant mass of the ZH system shows no shape variation. 
In this process, we only have initial state radiation and therefore significant changes in 
the shape are not expected for an observable like mzH■ Other observables, on the other 
hand, are highly sensitive to the choice of shower parameters. The distributions for the 
transverse momentum of the Higgs, p ^, but more importantly the transverse momentum 
of the ZH system, p!^ H , and that of the hardest jet, , which are trivially zero at parton- 
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Figure 8: Invariant ZH mass and distributions for gluon induced ZH production at ^fs = 14 
TeV. The left column shows the results obtained for the gg —> ZH case, with different starting 
scale for the shower: “wimpy” and “power” shower. The blue band shows the variation of the 
shower scale for “wimpy” shower in the range 0.5 pf < Qps < 2/xp, while the yellow bands show 
the uncertainty associated with a factor of two variation of the renormalisation and factorisation 
scales with respect to their central value. The right column shows the same results for the merged 
sample. The green curves in the left column correspond to the parton level results before passing 
them through Pythia 8. 


level, depend strongly on the shower parameters. We first notice that the shower produces 
a pff distribution harder than the parton-level one for all shower scale choices. This is 
related to the harder behaviour of the 2 —» 3 distributions discussed earlier. 

Another interesting observation to be made is related to the shape changes associated 
with the shower scale choice. The “power”-shower leads to consistently harder distribu¬ 
tions, while the “wimpy”-shower gives softer distributions. The different shower predictions 
start to diverge in a region correlated with the invariant mass of the ZH system, as this 
is the factorisation scale which is taken to be the starting scale of the “wimpy”-shower. 
The shower scale uncertainty bands become wider at larger pp values. This is more evi¬ 
dent in the second set of observables, p^ H and jP T , for which the non-merged predictions 
can vary by more than one order of magnitude between different shower scale options. 
At high transverse momentum, the shower uncertainty becomes more important than the 
intrinsic QCD one associated to the factorisation and renormalisation scale choice for the 
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Figure 9: Transverse momentum of the ZH system and hardest jet pt distributions for ZH 
production at y/s = 14 TeV. The setup is the same as in Fig. 8. 


hard process. We note here that despite the fact that the factorisation and renormalisation 
scale uncertainty is large, as evident from the yellow bands, it seems to mainly affect the 
normalisation of the curves. 

The advantage of the ME+PS procedure is then made obvious by noticing that the 
shower scale uncertainty is almost completely eliminated in the merged predictions. For 
all observables, the shower scale uncertainty bands remain well within the corresponding 
renormalisation and factorisation scale uncertainty ones, even at high transverse momen¬ 
tum. ME+PS predictions are therefore more accurate/precise and predictive than the 
parton shower alone as they include the exact 2 —3 matrix elements. These play an 
important role in the phase space regions populated by highly boosted objects which is 
often the case for LHC searches. 


3. Z(j) production in the 2HDM 

In the previous section we discussed gluon induced ZH production in the SM, employing 
the ME+PS merging method to improve the accuracy of the predictions for the differential 
distributions at the LHC. In this section, we will follow a similar approach for a beyond 
the SM scenario. The case we consider in this work is the 2HDM, as a minimal extension 
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of the SM [70]. The 2HDM extends the scalar sector of the SM by introducing a second 
SU(2)l doublet <1>2) which leads to five physical Higgs bosons, i.e. in the case of CP 
conservation, the light CP-even one, h°, a heavier CP-even one, H°, a CP-odd one, A 0 , 
and two charged Higgs bosons H^. Assuming no extra sources of CP-violation, seven 
input parameters fully characterise the model: 


tan /3, sin a, m/jO, m^o, m^o, m//±, m\ 2 , (3-1) 

with the convention 0 < (3 — a < it (with 0 < /3 < ir/2) fixing the sign of the Higgs 
coupling to the gauge bosons to be the same as in the SM. Depending on the structure of 
the Yukawa couplings, two main types of 2HDM setups can be considered: type / where 
all fermions couple to just one of the Higgs doublets and type II, where up-quarks (down) 
couple only to <f ?2 (3>i). 

Various theoretical requirements such as stability, perturbativity and unitarity, impose 
constraints on the 2HDM parameter space. At the same time, electroweak precision mea¬ 
surements and recent LHC Higgs physics results further constrain the parameter space, as 
discussed in more detail in [51]. Nevertheless, 2HDM scenarios that satisfy these constraints 
and yet have a significantly different phenomenology than the SM exist and have been stud¬ 
ied extensively in the literature. These scenarios arise in the “decoupling” limit [71], i.e., in 
the limit of cos(/3 — a) <C 1, which ensures that the masses of the additional Higgs bosons 
lie well above the light-Higgs one. Even in this case, significant shifts from the SM cou¬ 
plings are allowed, in particular in the tan/3 3> 1 limit, known as “delayed decoupling” [72]. 
Interestingly, the “decoupling” limit is not the only 2HDM realisation which is consistent 
with all parameter space constraints. Scenarios with light additional Higgs bosons are also 
allowed in the so-called “alignment limit” [71]. 

Some examples of viable 2HDM scenarios will be considered in this section to explore 
possible 2HDM signatures in Higgs production in association with a Z boson. Interesting 
features can arise in this process, not only because of possible deviations of the light Higgs 
couplings from their SM values, but most importantly because of the presence of the heavier 
states, H° and A 0 , which can lead to resonant production of Zcj) final states. Three neutral 
combinations of final states are possible: Zh°, ZH° and ZA°. These 2HDM processes have 
already been discussed in [29]. Similarly to ZH production in the SM, the production of 
the Zh° and ZH° final states can occur through Drell-Yan type diagrams, and in gluon- 
gluon fusion. The Drell-Yan like cross sections can be obtained straightforwardly by the 
appropriate rescaling of the SM cross-sections by the ratio of the g zz coupling over its 
SM counterpart, but the situation for the gluon fusion case is more involved. This can be 
inferred by considering the corresponding Feynman diagrams for the gluon fusion processes, 
shown in Figs. 0 and [ll]. The possibility of resonant production depends on the masses of 
A 0 and H°, while interesting interference patterns can arise due to relative sign of the A°cj)Z 
couplings. For completeness and to facilitate the discussion that follows, the dependance 
of the relevant Yukawa couplings on the 2HDM parameters is shown in Table |||, for type-/ 
and type-// setups, as rescalings of their SM counterparts. We note that the following 
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Coupling 

type-/ 

type-// 

9u 

cos a/ sin /3 

cos a/ sin /3 

-/i° 

9d 

cos a/ sin /3 

— since/ cos /3 

~H° 

9 U 

sin a/ sin/3 

sin a/ sin/3 

9d° 

since/ sin/3 

cos a/ cos /3 

9u° 

cot /3 

cot /3 

9d° 

— cot /3 

tan /3 


Table 3: Dependence of Yukawa couplings for up and down-type quarks on the 2HDM parameters 
for type-/ and type-// setups. The expressions in the table correspond to the ratio of the couplings 
over the corresponding SM value. 

couplings are also relevant for these process (valid for both type-/ and type-// setups): 

gy V = sm(/3 - a), gyy=cos{/3-a), gfy = 0, (3.2) 

with cjyy being the rescaling of the <f>VV coupling compared to the HVV one in the SM, 
while the A°(f>Z couplings are proportional to: 

gA°h° _ cos (jj _ and g f H ° — _ s in(/3 — a). (3-3) 
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Figure 10: Representative Feynman diagrams for ZH°/Zh° production in gluo 
2HDM. 
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Figure 11: Representative Feynman diagrams for ZA° production in gluon fusion in the 2HDM. 


We stress here that several studies have been presented in the literature in particular 
for the A°Z process, mostly in the context of the MSSM [73-77]. In the case of the 
MSSM, there are more constraints on the values of the Higgs couplings, while the 2HDM 
allows more freedom that can lead to more striking signals. A particularly interesting 
cosmologically motivated 2HDM scenario leading to a A 0 —> ZH° signature at the LHC is 
presented in [78], that finds very good prospects for discovery or exclusion even for the low- 
luminosity LHC. We also mention that various 2HDM scenarios allow significantly enhanced 
bottom Yukawas, and the Z(j) states can be produced mainly through bb annihilation. This 
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has been extensively discussed in the literature [76,79,80], and in relation with the subtleties 
of the treatment of the bottom quarks [81]. In this work, we will be focussing on the gluon 
fusion channel, presenting results for a series of 2HDM benchmarks. 


3.1 Calculation setup 


The calculation setup, regarding the reweighting and the ME+PS merging procedure, fol¬ 
lows closely that described in the previous section for the SM. In this section we discuss 
the details specific to the 2HDM implementation. The computation is performed within 
the MadGraph5_aMC@NLO framework [20], where MadLoop [47] is used for the com¬ 
putation of the one-loop amplitudes. The 2HDM@NLO model obtained from the pack¬ 
age NLOCT [82] is imported into MadGraph5_aMC@NLO for the computation of the 
2HDM amplitudes. The model is based on the FeynRules [83] and UFO [53] frameworks. 
More importantly for our computation, it includes all the necessary UV counterterms and 
R2 vertices for the MadLoop calculation. The model allows the computation of tree-level 
and one-loop amplitudes within a completely general 2HDM setup. The 2HDM parame¬ 
ters for the different benchmarks are imported into MadGraph5_aMC@NLO [20] using 
a parameter card, constructed with the help of an in-house modification of the public 
calculator 2HDMC [84], in the same way as in [51]. 

We stress again here that as described in [51] the 2HDM benchmarks to be used here 
have been constructed in agreement with all up-to-date parameter space constraints, which 
we have included through an interface of the public tools 2HDMC [84], HiggsBounds [85, 
86 ], SuperIso [87,88] and HiggsSignals [89,90] along with additional routines of our own. 

Three benchmarks will be employed to present the 2HDM results, with the correspond¬ 
ing parameters shown in Table [|. Benchmarks B1 and B2 have been constructed and used 
already for our study of Higgs pair production in the 2HDM [51], they correspond to B1 
and B4 in [51]. Benchmark B3 is a new one designed for this study. Here we briefly men¬ 
tion the main features of each benchmark. For completeness we also show the couplings 
relevant for gg —» Zcj) production in Table as rescalings of the SM couplings, similarly 
to Table 0 and Eqs. 3.2 -JO. 


• Benchmark Bl: A type-II 2HDM scenario with moderately heavy Higgs masses. 
Small tan/3 and cos(/3 — a) values ensure that the couplings of the light Higgs boson 
remain SM-like. The bottom Yukawa is slightly enhanced. This scenario allows a 
resonant production of both the light and Heavy Higgs with a Z boson through 
the decay of the pseudoscalar A 0 . The sign of the Zh°A° coupling determines the 
interference of the A°-mediated production with the SM-like diagrams. 

• Benchmark B2: A type-I 2HDM scenario with a relatively light heavy Higgs H° and 
a significantly heavier pseudoscalar A 0 . Both light-Higgs top and bottom Yukawas are 
enhanced by ~ 10%. The negative sign of m\ 2 protects the stability of the vacuum. 
This scenario also allows the resonant production of both the light and Heavy Higgs 
with a Z boson through the decay of the pseudoscalar A °. 
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tan /3 

a/n 

mu o 

m A o 

m H ± 

m 12 

B1 

1.75 

-0.1872 

300 

441 

442 

38300 

B2 

1.20 

-0.1760 

200 

500 

500 

-60000 

B3 

1.70 

-0.1757 

350 

250 

350 

12000 


Table 4: Parameter choices for the different 2HDM benchmarks used in our study. All masses are 
given in GeV. The lightest Higgs mass is fixed in all cases to m^o = 125 GeV. 



9h Q tt 

Sh°bb 

dH°tt 

9H°bb 

9A°tt 

9A°bb 

9A°Zh° 

9a°zh° 

9zzh° 

9zzh° 

B1 

0.958 

1.118 

-0.639 

1.677 

0.571 

1.75 

-0.069 

-0.998 

-0.0689 

0.998 

B2 

1.108 

1.108 

-0.684 

-0.684 

0.833 

-0.833 

0.141 

-0.990 

0.141 

0.990 

B3 

0.987 

1.034 

-0.608 

1.679 

0.588 

1.700 

-0.020 

- 1.000 

-0.020 

1.000 


Table 5: Normalised heavy-quark Yukawa couplings and Higgs Z couplings for the different 2HDM 
benchmarks defined in Table 3. Yukawa couplings are normalised to their SM counterparts as 
discussed in the text, while for the A°ZH° and A°ZH° couplings we show the proportionality 
constants of Eq. 3.3. 

• Benchmark B3: Another type-II 2HDM scenario with a reversed mass hierarchy 
between the heavy scalar H° and the pseudoscalar A 0 . The small tan/3 value allows 
us not to over-suppress the gA°tt coupling, while the c/A^bb is enhanced. Thanks to 
the inverted mass hierarchy m^o < m^o < o the resonant production of A 0 with 
a Z boson due to the heavy neutral Higgs decay becomes kinematically allowed. 

3.2 2HDM results 

In this section we present our results for the three 2HDM benchmarks introduced in the 
previous paragraph. We start by considering the total cross section for each process, which 
is shown in Table |6[ The heavy quark masses are again set to 173 and 4.75 GeV for top and 
bottom quarks, and the light Higgs mass to 125 GeV. The rest of the calculation details, 
such as the scale and PDF choices follow closely those of the SM calculation. We note here 
that where possible, we compared our results with the VH@NNLO version described in [29] 
and found very good agreement between the two implementations. 

Before moving to the discussion of some differential results, we first comment on the 
results in Table |6| First we notice that the cross-section for the Zh° process can be signif¬ 
icantly enhanced. To be more precise, benchmark B3 leads to a cross section nearly twice 
the SM prediction, benchmark B1 to a 60% enhancement, while B2 is gives a smaller ~20% 
increase. The main source of the increase in the cross-section is the presence of the resonant 
decay A 0 —>• ZhP, which is kinematically allowed in all three scenarios. The relative change 
in the ZhP cross section is strongly correlated with the mass of the pseudoscalar and the 
value of the A°Zh 0 coupling. We remind ourselves that this coupling is proportional to 
cos(/3 — a), i.e., it tends to zero in the alignment limit. For all scenarios considered here, 
its value remains small as seen in Table Consequently, it is not possible for this process 
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gg -+ Zh° 

gg -+ ZH° 
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B1 
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B3 

I 1 O +30% 
no _ 21% 

or q +30.1% 
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167 + 31 % 
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0.891 +33% 
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0.869 «« 

1 ooc +28% 
10 Z0 _ 21% 


Table 6: Cross sections (in fb) for gluon induced Z Higgs associated production at the LHC at 
•yds = 14 TeV for three 2HDM benchmarks. The uncertainties (in percent) refer to scale variations. 
No cuts are applied to final state particles and no Higgs or Z branching ratios are included. 


to receive extremely large contributions from the resonance. This is in contrast with what 
we have seen in light Higgs pair production where the resonant decay of the heavy Higgs 
can lead to an enhancement of up to a factor of 60 for the gg -+ h°h° cross section [51]. 

The most interesting feature of Table |6], is the potential size of the cross section for the 
ZH° process. We find that this can exceed 1 pb when the pseudoscalar A 0 is sufficiently 
heavy to allow the resonant decay into the heavy Higgs and a Z. This has been noticed and 
discussed recently in [78], as a signature for a cosmologically motivated 2HDM scenario. It 
is remarkable that even if the production threshold lies significantly higher, this process can 
lead to larger cross sections compared to the ZhP. This is possible as the relevant coupling, 
ZH°A°, as shown in Table ||, is not suppressed by the “SM-like” light Higgs constraints. 
Despite the fact that the prospects for discovery depend strongly on the resulting decay 
products of the heavy Higgs, it is worth noting that even in the scenarios where H° decays 
predominantly into bb, the current experimental searches for ZH set a cut on the invariant 
mass of the bb pair close to the light Higgs mass and would therefore miss this signal. 
Finally, we note that the ZA° production cross section remains very small in the scenarios 
where the A 0 is heavier than H°, but can reach the picobarn level in a scenario such as 
benchmark B3, as a result of the inverted mass hierarchy. 

Further interesting information on these processes can be extracted from the differential 
distributions. For brevity we present only those for the invariant mass of the system and the 
transverse momentum of the Higgs, but our setup is fully differential and any distribution 
can be plotted. We show these in Fig. 12, for the cases in which the cross section is not 
negligible. The results shown here are obtained with merged samples of 0 and 1-jet matched 
to Pythia 8 for parton shower, in the same setup as that described in Section 2 for the 
SM. 


For the Zh° final state we also show the SM prediction for comparison. Resonance 
peaks arise in all scenarios for Zh°, each time located at the mass of the pseudoscalar 
A 0 . The sharpness of the peak varies with the mass of H°, as heavier have larger 
widths going from 0.01 GeV for B3, to 7 GeV for B1 and 35 GeV in B2. We also notice 
various interesting interference patterns, clearly visible for benchmarks B1 and B2. The 
H°-mediated diagram interferes with the SM-like amplitude, with the interference switching 
sign at Comparing scenarios B1 and B2, we see that the Zh°A° couplings have 

opposite signs and therefore in one case the dip appears right before the resonance peak, 
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Figure 12: Differential distributions for gluon induced Zh°, ZH° and ZA° production at y/s = 14 
TeV for the three 2HDM benchmarks and comparison with the SM. Left: Invariant mass of the Z(f> 
system. Right: Transverse momentum of the Higgs. 


and in the other right after. More subtle features are also visible in the plots away from the 
resonance peaks. These features can always be traced back to the 2HDM parameters and 
the value of the relevant couplings as shown in Table |5|. One such example is the fact that 
the B2 mzh 0 curve lies a bit lower than the SM one in the region below 350 GeV, which 
can be linked to the enhanced top Yukawa leading to a bigger box contribution. The box 
is in turn interfering destructively with the triangle leading to a smaller total amplitude 
for the gg —>• Zh° process. 
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For the ZH° process, only the two benchmarks that give measurable cross sections are 
shown. The plots shown for this process are dominated by the resonant decay of A 0 . This 
is more obvious in the B1 curve as the resonant peak is closer to the threshold. Scenario B2 
receives some non-negligible off-peak contributions from the Z triangle and box diagrams, 
which in this case interfere constructively, as the H° top Yukawa sign is flipped. For B1 
both the top Yukawa and ZZH° couplings signs are flipped, therefore the interference 
between triangle and box is destructive, and the result in the tails away from the resonant 
peaks, is suppressed compared to B2. 

The situation is less complicated for ZA° for which in B3 a resonance very close to 
the mz + itia 0 threshold dominates the plots, while the cross sections for B1 and B2 
are extremely suppressed as no resonant decay is kinematically allowed. Moreover the 
production of a rather heavy ZA° pair probes the gluon luminosity at large partonic x 
values and is therefore suppressed. 


4. Conclusions 

Investigating the nature of the Higgs boson discovered at the LHC is a challenging task. 
While the results of the measurements undertaken so far show that the 125 GeV scalar 
agrees well with the SM prediction, there is still room for deviations from the SM and 
possibly an extended Higgs sector to be discovered at the LHC. The exploration of various 
Higgs production processes is of vital importance to exclude or confirm a non-minimal 
Higgs sector. 

An important process yet to be precisely measured at the LHC is the associated pro¬ 
duction of a Higgs with a Z boson. In addition to the Drell-Yan type contributions, this 
process acquires a gluon fusion component at NNLO, which proves to be of particular im¬ 
portance in the boosted regime. In this work, we have reviewed the main features of the 
gg —>• ZH process, both at the matrix-element and cross-section level. We have examined 
the behaviour and the relative importance of the 2 —> 2 and 2 —> 3 matrix elements for the 
gluon induced component. We have found that in the high px regions the 2—^3 matrix 
elements behave in a different way from the 2 —> 2 ones and therefore have to be taken 
into account to provide accurate predictions for the differential distributions. To achieve 
this, we have combined the two in a consistent way, by merging different jet multiplicity 
samples and matching them to a parton shower. 

Our results have been obtained within the MadGraph5_aMC@NLO framework with 
the help of Pythia 8 for the parton shower. The ME+PS approach provides a more ac¬ 
curate description of the process compared to the parton shower alone. In particular, it 
significantly reduces the uncertainty associated with the shower scale choice. For observ¬ 
ables such as the transverse momentum of radiated jets in the hard region, the prediction of 
the parton shower alone can be misleading as here the results are extremely sensitive to the 
shower parameters. We find that in the merged predictions this sensitivity is almost com¬ 
pletely eliminated, with the shower uncertainty remaining well within the intrinsic QCD 
uncertainty due to the renormalisation and factorisation scale variations. 
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The reduction of the uncertainties associated with the SM prediction and especially the 
accurate description of differential distributions is crucial for searches for beyond the SM 
scenarios. One scenario that the LHC aims to explore is the 2HDM. In this paper, we have 
also provided predictions for the gluon fusion component of the Zcp associated production 
in the 2HDM. Following the same setup as in the SM, we have presented our predictions 
for three representative 2HDM benchmarks. We have considered all three neutral Higgs 
bosons, presenting results for the cross sections and the differential distributions. 

In the production of the light Higgs in association with a Z, large enhancements can 
be achieved compared to the SM prediction if the resonant decay of the pseudoscalar A 0 
is kinematically allowed. Moreover, interference patterns arise between the additional dia¬ 
grams and the SM-like ones, leading to interesting features in the differential distributions. 
The resonant production of a H°Z pair also becomes important as the H°ZA° coupling 
is not suppressed, leading to large cross sections for gg —>• ZH° if the pseudoscalar A 0 
is heavier then H°. Finally in scenarios where the pseudoscalar A 0 is lighter than the 
heavy Higgs, gg —>• H° —>• ZA° production is allowed and leads to large cross sections in 
still-to-be-excluded scenarios. 


Acknowledgements 

We would like to thank David Lopez-Val for his assistance with the 2HDM benchmarks, 
Stefan Liebler for a private version of vh@nnlo with the 2HDM implementation and Am- 
bresh Shivaji and Olivier Mattelaer for useful discussions. This work has been performed in 
the framework of the ERC Grant No. 291377 “LHCTheory” and has been supported in part 
by the European Union as part of the FP7 Marie Curie Initial Training Network MCnetITN 
(PITN-GA-2012-315877) and by the National Fund for Scientific Research (F.R.S.-FNRS 
Belgium) under a FRIA grant. 


References 


[1] ATLAS Collaboration, G. Aad et al., Observation of a new particle in the search for the 
Standard Model Higgs boson with the ATLAS detector at the LHC, Phys.Lett. B716 (2012) 
1-29, JarXiv: 1207.7214J . 

[2] CMS Collaboration, S. Chatrchyan et al., Observation of a new boson at a mass of 125 GeV 
with the CMS experiment at the LHC, Phys.Lett. B716 (2012) 30-61, arXiv: 1207.7235 |. 

[3] S. Weinberg, Implications of Dynamical Symmetry Breaking, Phys.Rev. D13 (1976) 974-996. 

[4] F. Englert and R. Brout, Broken Symmetry and the Mass of Gauge Vector Mesons, 

Phys.Rev.Lett. 13 (1964) 321-323. 

[5] P. W. Higgs, Broken Symmetries and the Masses of Gauge Bosons, Phys.Rev.Lett. 13 (1964) 
508-509. 

[6] TEVNPH (Tevatron New Phenomena and Higgs Working Group), CDF 
Collaboration, DO Collaboration, Combined CDF and DO Search for Standard Model Higgs 
Boson Production with up to 10.0 fb -1 of Data, 


arXiv:1203.3774. 


- 23 - 





[7] ATLAS Collaboration, G. Aad et al., Search for the bb decay of the Standard Model Higgs 
boson in associated ( W/Z)H production with the ATLAS detector, arXiv: 1409.6212. 


[ 8 ] 

[9] 


[ 10 ] 

[ 11 ] 

[ 12 ] 

[13] 

[14] 

[15] 

[16] 
[17] 


[18] 

[19] 


[ 20 ] 


[ 21 ] 


[ 22 ] 

[23] 


CMS Collaboration, S. Chatrchyan et al., Search for the standard model Higgs boson 
produced in association with a W or a Z boson and decaying to bottom quarks, Phys.Rev. 
D89 (2014), no. 1 012003, 


arXiv:1310.3687 


Search for associated production of the Higgs boson in the WH —> WWW^ —» lulvlu and 
ZH —> ZWW^*l — > lllvlv channels with the ATLAS detector at the LHC, Tech. Rep. 
ATLAS-CONF-2013-075, CERN, Geneva, Jul, 2013. 


CMS Collaboration, VH with H—> WW—> ivtv and V—¥jj, Tech. Rep. CMS-PAS-HIG-13-017, 
CERN, Geneva, 2013. 


ATLAS Collaboration, Search for invisible decays of a Higgs boson produced in association 
with a Z boson in ATLAS, Tech. Rep. ATLAS-CONF-2013-011, CERN, Geneva, Mar, 2013. 


CMS Collaboration, Search for invisible Higgs produced in association with a Z boson, Tech. 
Rep. CMS-PAS-HIG-13-018, CERN, Geneva, 2013. 

R. Hamberg, W. van Neerven, and T. Matsuura, A Complete calculation of the order or s 
correction to the Drell-Yan K factor, Nucl.Phys. B359 (1991) 343-405. 


O. Brein, A. Djouadi, and R. Harlander, NNLO QCD corrections to the Higgs-strahlung 
processes at hadron colliders, Phys.Lett. B579 (2004) 149-156, [hep-ph/0307206]. 


O. Brein, R. Harlander, M. Wiesemann, and T. Zirke, Top-Quark Mediated Effects in 


Hadronic Higgs-Strahlung, Eur.Phys.J. C72 (2012) 1868, [arXiv: 1111.0761]. 


O. Brein, R. V. Harlander, and T. J. Zirke, vh@nnlo - Higgs Strahlung at hadron colliders, 


Comput.Phys.Commun. 184 (2013) 998-1003, [arXiv : 1210.5347]. 


G. Ferrera, M. Grazzini, and F. Tramontano, Associated ZH production at hadron colliders: 
the fully differential NNLO QCD calculation, Phys.Lett. B740 (2015) 51-55, 

[ arXiv: 1407.4747 ]. 


M. Ciccolini, S. Dittmaier, and M. Kramer, Electroweak radiative corrections to associated 
WH and ZH production at hadron colliders, Phys.Rev. D68 (2003) 073003, [hep-ph/0306234]. 

A. Denner, S. Dittmaier, S. Kallweit, and A. Mck, HAWK 2.0: A Monte Carlo program for 
Higgs production in vector-boson fusion and Higgs strahlung at hadron colliders, 

arXiv: 1412.539C| . 

J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, et al., The automated computation 
of tree-level and next-to-leading order differential cross sections, and their matching to parton 
shower simulations, JHEP 1407 (2014) 079, arXiv: 1405.0301]. 

S. Alioli, P. Nason, C. Oleari, and E. Re, A general framework for implementing NLO 
calculations in shower Monte Carlo programs: the POWHEG BOX, JHEP 1006 (2010) 043, 

| arXiv: 1002.2581 ], 


arXiv:1405.0301 


hep-ph/0306234 


T. Sjostrand, S. Mrenna, and P. Z. Skands, A Brief Introduction to PYTHIA 8.1, 


Comput.Phys.Commun. 178 (2008) 852-867, [arXiv:0710.3820|. 


M. Bahr, S. Gieseke, M. Gigg, D. Grellscheid, K. Hamilton, et al., Herwig++ Physics and 


Manual, Eur.Phys.J. C58 (2008) 639-707, [arXiv : 0803.0883]. 


- 24 
















[24] T. Gleisberg, S. Hoeclie, F. Krauss, M. Schonherr, S. Schumann, et al., Event generation with 
SHERPA 1.1, JHEP 0902 (2009) 007, jarXiv: 0811.462^ . 

[25] D. A. Dicus and C. Kao, Higgs Boson - Z° Production From Gluon Fusion, Phys.Rev. D38 
(1988) 1008. 

[26] B. A. Kniehl, Associated Production of Higgs and Z Bosons From Gluon Fusion in Hadron 
Collisions, Phys.Rev. D42 (1990) 2253-2258. 

[27] L. Altenkamp, S. Dittmaier, R. V. Harlander, H. Rzehak, and T. J. Zirke, Gluon-induced 


Higgs-strahlung at next-to-leading order QCD, JHEP 1302 (2013) 078, [arXiv: 1211.5015]. 

[28] R. V. Harlander, A. Kulesza, V. Tlieeuwes, and T. Zirke, Soft gluon resummation for 
gluon-induced Higgs Strahlung, JHEP 1411 (2014) 082, [arXiv: 1410.0217|. 

[29] R. V. Harlander, S. Liebler, and T. Zirke, Higgs Strahlung at the Large Hadron Collider in 
the 2-Higgs-Doublet Model, JHEP 1402 (2014) 023, [arXiv: 1307.8122]. 

[30] CMS Collaboration, Search for a pseudoscalar boson A decaying into a Z and an h boson in 
the llbb final state, Tech. Rep. CMS-PAS-HIG-14-011, CERN, Geneva, 2014. 

[31] CMS Collaboration, 2HDM scenario, H to h.h and A to Zh, Tech. Rep. 
CMS-PAS-HIG-13-025, CERN, Geneva, 2013. 

[32] C. Englert, M. McCullough, and M. Spannowsky, Gluon-initiated associated production boosts 
Higgs physics, Phys.Rev. D89 (2014), no. 1 013013, [arXiv: 1310.4828|. 

[33] J. M. Butterworth, A. R. Davison, M. Rubin, and G. P. Salarn, Jet substructure as a new 
Higgs search channel at the LHC, Phys.Rev.Lett. 100 (2008) 242001, | arXiv : 0802. 24~7ct] . 

[34] J. Alwall, Q. Li, and F. Maltoni, Matched predictions for Higgs production via heavy-quark 
loops in the SM and beyond, Phys.Rev. D85 (2012) 014031, [arXiv: 1110.1728]. 

[35] Q. Li, Q.-S. Yan, and X. Zhao, Higgs Pair Production: Improved Description by Matrix 
Element Matching, Phys.Rev. D89 (2014) 033015, [arXiv: 1312.383C]. 

[36] P. Maierhofer and A. Papaefstathiou, Higgs Boson pair production merged to one jet, JHEP 
1403 (2014) 126, farXiv: 1401.0007| ]. 

[37] S. Catani, F. Krauss, R. Kuhn, and B. Webber, QCD matrix elements + parton showers, 
JHEP 0111 (2001) 063, flh7p- ph/010923l |. 

[38] F. Krauss, Matrix elements and parton showers in hadronic interactions, JHEP 0208 (2002) 
015, fhep - ph/0205283 |. 

[39] L. Lonnblad, Correcting the color dipole cascade model with fixed order matrix elements, 
JHEP 0205 (2002) 046, ]^p- ph/0112284 |. 

[40] N. Lavesson and L. Lonnblad, W+jets matrix elements and the dipole cascade, JHEP 0507 
(2005) 054, fhep-p h/0503293 |. 

[41] S. Hoeche, F. Krauss, S. Schumann, and F. Siegert, QCD matrix elements and truncated 
showers, JHEP 0905 (2009) 053, [)arXiv:0903. 12~19 |. 

[42] K. Hamilton, P. Richardson, and J. Tully, A Modified CKKW matrix element merging 
approach to angular-ordered parton showers, JHEP 0911 (2009) 038, ||arXiv: 0905.3072 1. 

[43] M. L. Mangano, M. Moretti, F. Piccinini, and M. Treccani, Matching matrix elements and 
shower evolution for top-quark production in hadronic collisions, JHEP 0701 (2007) 013, 

| hep- p h/0611129 |. 


- 25 - 


























[44] R. Frederix and S. Frixione, Merging meets matching in MC@NLO, JHEP 1212 (2012) 061, 


arXiv:1209.6215 . 


[45] L. Lonnblad and S. Prestel, Merging Multi-leg NLO Matrix Elements with Parton Showers, 
JHEP 1303 (2013) 166, flarXiv: 1211.7278 |. 

[46] P. Agrawal and A. Shivaji, Gluon Fusion Contribution to VHj Production at Hadron 
Colliders, Phys.Lett. B741 (2014) 111-116, [arXiv : 1409.80591. 

[47] V. Hirsclii et al., Automation of one-loop QCD corrections, JHEP 05 (2011) 044, 

| arXiv: 1103.0621 ]. 

[48] G. Ossola, C. G. Papadopoulos, and R. Pittau, Reducing full one-loop amplitudes to scalar 
integrals at the integrand level, Nucl.Phys. B763 (2007) 147-169, |tiep-ph/0609007]. 

[49] G. Ossola, C. G. Papadopoulos, and R. Pittau, CutTools: A Program implementing the OPP 
reduction method to compute one-loop amplitudes, JHEP 0803 (2008) 042, 

| arXiv:0711.3596 ]. 

[50] V. Hirsclii and O. Mattelaer. in preparation. 

[51] B. Hespel, D. Lopez-Val, and E. Vryonidou, Higgs pair production via gluon fusion in the 
Two-Higgs-Doublet Model , JHEP 1409 (2014) 124, flarXi v : 1407.0281 1. 

[52] F. Maltoni, E. Vryonidou, and M. Zaro, Top-quark mass effects in double and triple Higgs 
production in gluon-gluon fusion at NLO, JHEP 1411 (2014) 079, [jarXiv : 1408.6542 1 . 

[53] C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer, et al., UFO - The Universal 
FeynRules Output, Comput.Phys.Commun. 183 (2012) 1201 1214, [arXiv:1108.2040 . 

[54] P. de Aquino, W. Link, F. Maltoni, O. Mattelaer, and T. Stelzer, ALOHA: Automatic 
Libraries Of Helicity Amplitudes for Feynman Diagram Computations, 
Comput.Phys.Commun. 183 (2012) 2254-2263, [arXiv:1108.2041 . 

[55] A. Denner, S. Dittmaier, M. Roth, and D. Wackeroth, Predictions for all processes e+ e- 4 
fermions + gamma, Nucl.Phys. B560 (1999) 33-65, hep-ph/990447^| . 

[56] A. Denner, S. Dittmaier, M. Roth, and L. Wieders, Electroweak corrections to 

charged-current e+ e- 4 fermion processes: Technical details and further results, Nucl.Phys. 
B724 (2005) 247-294, [^p- ph/0505042 |. 

[57] C. Anastasiou, S. Buehler, F. Herzog, and A. Lazopoulos, Total cross-section for Higgs boson 
hadroproduction with anomalous Standard Model interactions, JHEP 1112 (2011) 058, 

| arXiv: 1107.0683)1 . 

[58] A. Martin, W. Stirling, R. Thorne, and G. Watt, Parton distributions for the LHC, 
Eur.Phys.J. C63 (2009) 189-285, flarXiv: 0901.0002 ], 

[59] R. Frederix, S. Frixione, V. Hirsclii, F. Maltoni, R. Pittau, et al., Four-lepton production at 
hadron colliders: aMC@NLO predictions with theoretical uncertainties, JHEP 1202 (2012) 
099, farXiv: 1110.4738|] . 

[60] S. Frixione and B. R. Webber, Matching NLO QCD computations and parton shower 
simulations, JHEP 06 (2002) 029, [tiep-ph/0204244|. 

[61] F. Maltoni, K. Paul, T. Stelzer, and S. Willenbrock, Associated production of Higgs and 
single top at hadron colliders, Phys.Rev. D64 (2001) 094023, [tiep-ph/0106293|. 


- 26 - 



























[62] S. Biswas, E. Gabrielli, and B. Mele, Single top and Higgs associated production as a probe of 
the Htt coupling sign at the LHC, JHEP 1301 (2013) 088, | arXiv: 1211. 049S ]. 

[63] M. Farina, C. Grojean, F. Maltoni, E. Salvioni, and A. Thamm, Lifting degeneracies in Higgs 
couplings using single top production in association with a Higgs boson, JHEP 1305 (2013) 
022, JarXiv: 1211,3736|] . 

[64] F. Maltoni, K. Mawatari, and M. Zaro, Higgs characterisation via vector-boson fusion and 
associated production: NLO and parton-shower effects, Eur.Phys.J. C74 (2014), no. 1 2710, 
[ arXiv: 1311.1829 ]. 

[65] W. Stirling and E. Vryonidou, Electroweak gauge boson polarisation at the LHC, JHEP 1207 
(2012) 124, JarXiv: 1204.6427 ]. 

[66] Z. Bern, G. Diana, L. Dixon, F. Febres Cordero, D. Forde, et al., Left-Handed W Bosons at 
the LHC, Phys.Rev. D84 (2011) 034008, [arXiv: 1103.5445], 

[67] T. Sjostrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, et al., An Introduction to 


PYTHIA 8.2, 


arXiv:1410.3012. 


[68] J. Alwall, S. de Visscher, and F. Maltoni, QCD radiation in the production of heavy colored 
particles at the LHC, JHEP 0902 (2009) 017, [arXiv:0810.5350|. 

[69] P. de Aquino, F. Maltoni, K. Mawatari, and B. Oexl, Light Gravitino Production in 
Association with Gluinos at the LHC, JHEP 1210 (2012) 008, [arXiv: 1206.7098]. 

[70] G. Branco, P. Ferreira, L. Lavoura, M. Rebelo, M. Sher, et al., Theory and phenomenology of 
two-Higgs-doublet models, Phys.Rept. 516 (2012) 1-102, [arXiv: 1106.0034|. 

[71] J. F. Gunion and H. E. Haber, The CP conserving two Higgs doublet model: The Approach to 
the decoupling limit, Phys.Rev. D67 (2003) 075019, [iep-ph/020701C]. 

[72] H. E. Haber, M. J. Herrero, H. E. Logan, S. Penaranda, S. Rigolin, et al., SUSY QCD 
corrections to the MSSM hO bb vertex in the decoupling limit, Phys.Rev. D63 (2001) 055004, 


[hep- ph/0007006 


[73] C. Kao, Production of a pseudoscalar Higgs with a Z boson from gluon fusion, Phys.Rev. D46 
(1992) 4907-4913. 

[74] J. Yin, W.-G. Ma, R.-Y. Zhang, and H.-S. Hou, A0Z0 associated production at the CERN 
large hadron collider in the minimal supersymmetric standard model, Phys.Rev. D66 (2002) 
095008. 

[75] C. Kao, G. Lovelace, and L. H. Orr, Detecting a Higgs pseudoscalar with a Z boson at the 
LHC, Phys.Lett. B567 (2003) 259-264, [iep-ph/0305028|. 

[76] L. L. Yang, C. S. Li, J. J. Liu, and L. G. Jin, Production of scalar Higgs bosons associated 
with Z0 boson at the CERN LHC in the MSSM, J.Phys. G30 (2004) 1821-1835, 


[hep- ph/0312179 


[77] B. A. Kniehl and C. P. Palisoc, Associated production of Z and neutral Higgs bosons at the 


CERN Large Hadron Collider, Phys.Rev. D85 (2012) 075027, arXiv: 1112.1575 


[78] G. Dorsch, S. Huber, K. Mirnasu, and J. No, Echoes of the Electroweak Phase Transition: 
Discovering a second Higgs doublet through Aq —> ZHq, Phys.Rev.Lett. 113 (2014), no. 21 
211802, jarXiv: 1405.5537 ], 


- 27 - 



















[79] C. Kao and S. Sacliithanandam, Detecting a Higgs pseudoscalar with a Z boson produced in 
bottom quark fusion, Phys.Lett. B620 (2005) 80-87, | h.ep-ph/0411331 |. 

[80] Q. Li, C. S. Li, J. J. Liu, L. G. Jin, and C.-P. Yuan, Next-to-leading order QCD predictions 
for AO ZO associated production at the CERN large hadron collider, Phys.Rev. D72 (2005) 
034032, ihep-ph/0501070 ]. 

[81] R. Harlander, M. Kramer, and M. Schumacher, Bottom-quark associated Higgs-boson 
production: reconciling the four- and five-flavour scheme approach, |arXiv : 1112.3478. 


[82] C. Degrande, Automatic evaluation of UV and R2 terms for beyond the Standard Model 
Lagrangians: a proof-of-principle, 


arXiv:1406.3030 


[83] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, FeynRules 2.0 - A 
complete toolbox for tree-level phenomenology, Comput.Phys.Commun. 185 (2014) 2250-2300, 

[ arXiv: 1310.1921 ]. 

[84] D. Eriksson, J. Rathsman, and O. Stal, 2HDMC: Two-Higgs-Doublet Model Calculator 
Physics and Manual, Comput.Phys.Commun. 181 (2010) 189-205, [arXiv:0902.08511. 

[85] P. Beclitle, O. Brein, S. Heinemeyer, G. Weiglein, and K. E. Williams, HiggsBounds: 
Confronting Arbitrary Higgs Sectors with Exclusion Bounds from LEP and the Tevatron, 
Comput.Phys.Commun. 181 (2010) 138-167, | arXiv:08 1 1.4169 ], 

[86] P. Beclitle, O. Brein, S. Heinemeyer, G. Weiglein, and K. E. Williams, HiggsBounds 2.0.0: 
Confronting Neutral and Charged Higgs Sector Predictions with Exclusion Bounds from LEP 
and the Tevatron, Comput.Phys.Commun. 182 (2011) 2605-2631, [arXiv: 1102.18981. 

[87] F. Mahmoudi, Superlso: A Program for calculating the isospin asymmetry of B — » K *7 in 
the MSSM, Comput.Phys.Commun. 178 (2008) 745-754, | arXiv:0710.2067 ], 

[88] F. Mahmoudi, Superlso v2.3: A Program for calculating flavor physics observables in 
Supersymmetry, Comput.Phys.Commun. 180 (2009) 1579-1613, [ arXiv:0808.3144 ]. 

[89] P. Beclitle, S. Heinemeyer, O. Stal, T. Stefaniak, and G. Weiglein, HiggsSignals: 
Confronting arbitrary Higgs sectors with measurements at the Tevatron and the LHC, 
Eur.Phys.J. C74 (2014) 2711, flarXiv: 1305.1933j] . 

[90] O. Stal and T. Stefaniak, Constraining extended Higgs sectors with HiggsSignals, PoS 

EPS-HEP2013 (2013) 314, farXiv: 1310.4039|1. 


- 28 - 


















